
nl=require "luaNumeric"
local unity=require "luaUnity"
require "Q2DViewerEx"
local mpi = require "luaMPI"


--===========================================================================--
--初始化参数
local SBRSolver={

	path = ' ';                           --主目录:                c:/almond/
	projectPath = '';                  -- 总工程文件目录:  c:/almond/almond.gap
	resultsPath = 'MonoRCS';     -- 计算结果目录:      c:/almond/almond-RCS/
	gapsPath = 'gaps';               -- 计算结果目录:      c:/almond/almond-GAPs/
	geometryFilePath = '';          -- 导入模型目录:      c:/almond.stp
	fileName = 'rcs';                  --文件名:                c:/almond-details
	name = '';                           --文件名:                almond
	rcsProjectStem = '';              --
	name = 'RCS';
	
	bMonoRCS = true;
	
	bReserveTraces = false;

	bFilter = false;
	bInnerFilter = true;

	minFreq = 1;
	maxFreq = 5;
	numFreq = 6;
	bFreqList = false;
	freqList = '(1,10,10)';
	
	numThreads = 20;
	minMesh = 3;

	minTheta= 0;
	maxTheta = 180;
	numTheta= 181;
	bThetaList = false;
	thetaList = '(-90,90,91)';
	
	minPhi= 0;
	maxPhi = 0;
	numPhi= 1;
	bPhiList = false;
	phiList = '(0,90,3)';

	polarHH = true;
	polarVV = false;
	
	range_min = -10;
	range_max = 10;
	auto_range = false;
	
	scale_type = 1;
	
	data_post_fun = function(d) 
		if d == nil then d=1e-20; end
		local d = nl.abs(d)*nl.abs(d)*4.0*pi;
		if d < 1e-20 then d=1e-20; end
		return 10*nl.log10(d);
	end;
	
	results = {};
	
	set_settings_str = Gallop.set_settings_str;
	get_settings_str = Gallop.get_settings_str;
};
--===========================================================================--


--===========================================================================--
--创建表和创建元表
function SBRSolver:new (o)
	o = o or {};
	self.__index = self;	
	setmetatable(o,self)
	o:regist_funcs_to_gallop();
	o.RCS = nl.new_eigen_ctensor4(1,1,1,1);
	o.init_strs = o:get_settings_str();
	o:init();
	return o;
end
--===========================================================================--

--===========================================================================--
--初始化函数
function SBRSolver:init ()
	self.funs = {};
	self.funs.abs  = function(d) 
		if d == nil then d=1e-20; end
		local d = nl.abs(d)*nl.abs(d)*4.0*pi;
		if d < 1e-20 then d=1e-20; end
		return d;
	end;
	self.funs.phs = function(input) return nl.arg(input)*180/pi; end;
	self.funs.fields   = function(d) 
		if d == nil then d=1e-20; end
		local d = nl.abs(d);
		if d < 1e-20 then d=1e-20; end
		return d;
	end;
	self.funs.log  = function(d) 
		if d == nil then d=1e-20; end
		local d = nl.abs(d)*nl.abs(d)*4.0*pi;
		if d < 1e-20 then d=1e-20; end
		return 10*nl.log10(d);
	end;
	self.fun_names = {};
	self.fun_names[1] = 'log';
	self.fun_names[2] = 'abs';
	self.fun_names[3] = 'fields';
	self.fun_names[4] = 'phs';
	
	self.ylable_names = {};
	self.ylable_names[1] = 'RCS [dBsm]';
	self.ylable_names[2] = 'RCS [sm]';
	self.ylable_names[3] = 'Farfield [V/m]';
	self.ylable_names[4] = 'Farfield [arc.]';
end
--===========================================================================--

--===========================================================================--
--注册全局函数
function SBRSolver:regist_funcs_to_gallop()
	Gallop.init_funcs['SBRApp'] = function()
		local str = self.init_strs;
		assert( type('str')=='string', 'SBRApp Plugin init string error' );
		if str ~= '' then
			Gallop.SBRApp:set_settings_str(str)
		end
		Gallop.SBRApp:init_parameters();
		Gallop.SBRApp.RCS = nil;
		Gallop.SBRApp.shape = nil;
		Gallop.SBRApp.fields = {};	
		Gallop.SBRApp:update_post_fun();
		Gallop.SBRApp:plot();
		Gallop.SBRApp:update_shape();
		Gallop.SBRApp:update_border_to_ui();
		--Gallop.SBRApp:plot2D();
		Gallop.SBRApp:update_parameters_to_ui();
		Gallop.SBRApp:update_freqlist_to_ui();
		Gallop.SBRApp:plot2D();
		Gallop.SBRApp:update_view_rcs_type_to_ui();
		
		--Gallop.SBRApp:update_freq_index_to_ui();
	end;
	
	Gallop.load_funcs['SBRApp'] = function()
		local str = Geometry.get_setting( 'SBRApp.Settings' )
		if str ~= '' then
			Gallop.SBRApp:set_settings_str(str)
		end
		Gallop.SBRApp:init_parameters();
		Gallop.SBRApp:update_parameters_to_ui();
		Gallop.SBRApp:init_parameters()
		Gallop.SBRApp:update_post_fun();
		Gallop.SBRApp:plot();
		Gallop.SBRApp.shape = Geometry.object('SBR');
		Gallop.SBRApp:update_border_to_ui();
		Gallop.SBRApp:update_freqlist_to_ui();	
		--print( Gallop.SBRApp.freq_idx )
		Gallop.SBRApp:update_freq_index_to_ui();
		Gallop.SBRApp:plot2D();
		Gallop.SBRApp:update_view_rcs_type_to_ui();
	end;
	
	Gallop.presave_funcs['SBRApp'] = function()
		Gallop.SBRApp:update_parameters_from_ui()
		local str = Gallop.SBRApp:get_settings_str();
		Geometry.set_setting( 'SBRApp.Settings', str );
	end;
end
--===========================================================================--


--===========================================================================--
-- 保存和读取函数
function SBRSolver:set_settings_str(str)
	Gallop.set_settings_str(self,str);
	self.results = undump_t(self.strResults);
	if self.strRCS then
		self.RCS = undump_t(self.strRCS);
	end
	self.strRCS = nil;
	self.strResults = nil;
end

function SBRSolver:get_settings_str()
	self.strResults = dump_t(self.results,'');
	if self.RCS then 
		self.strRCS= nl.serialize(self.RCS);
	end
	local str = Gallop.get_settings_str(self);
	self.strResults = nil;
	self.strRCS = nil;
	return str;
end
--===========================================================================--

--===========================================================================--
--ui交互函数
function SBRSolver:update_parameters_from_ui()
end

function SBRSolver:update_parameters_to_ui()
end
function SBRSolver:update_border_to_ui()
end
function SBRSolver:update_results_path_to_ui()
end

function SBRSolver:update_freqlist_to_ui()
end

function SBRSolver:update_progress_bar(per)
end

function SBRSolver:update_freq_index_to_ui()
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:get_freq_list(list)
	list:clear();
	for k,v in pairs( self.probe_freq ) do 
		for p,pp in pairs( self.polar_angles ) do 
			list:add( '频率：'..v..' GHz ' .. '极化:' .. self.polars[p] );
		end
	end
	--print(list)
end
--===========================================================================--

--=====================================================================--
--更新画图函数
function SBRSolver:update_range()
	local fields = self.fields[self.freq_idx];
	if fields and self.auto_range then
		local range = { fields:value_range() };
		self.range_min = range[1];
		self.range_max = range[2];
	end
	self:update_range_to_ui();
	self:update_range_to_view();
end
--=====================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:change_selected_freq(idx)
	assert( idx <= #self.probe_freq * #self.polar_angles and idx >= 0 , '频率超出范围')
	if self.freq_idx ~= idx+1 then
		self.freq_idx = idx+1;
		self:update_freq_index_to_ui();
		self:update_fields_to_view();
		self:update_range()
		self:update_view();
	end
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:change_has_filter(set)
	--print( self.bFilter )
	self.bFilter = set;
end

function SBRSolver:change_is_inner_filter(set)
	--print( self.bInnerFilter )
	self.bInnerFilter = set;
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:change_scale_type(set)
	--print( self.bFilter )
	assert( set <=4 and set >=1, '类型设置错误' )
	if self.scale_type ~= set then 
		self.scale_type = set;
		self.data_post_fun = self.funs[ self.fun_names[set] ];
		self:plot();
		self:update_view_rcs_type_to_ui()
	end
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:update_post_fun()
	assert( self.scale_type <=4 and self.scale_type >=1, '类型设置错误' )
	print('fun name:',self.fun_names[self.scale_type])
	self.data_post_fun = self.funs[ self.fun_names[self.scale_type] ];
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:get_border()
	if self.shape then 
		--print(self.shape)
		self.bbox= self.shape:bbox();
		--print(self.bbox)
		self.bbox:extend(self.bbox:size()/2)
		return self.bbox:xmin(),self.bbox:ymin(),self.bbox:zmin(),self.bbox:xmax(),self.bbox:ymax(),self.bbox:zmax();
	end
	return -1,-1,-1,1,1,1;
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:change_view_range(min,max,bset)
	self.range_min,self.range_max,self.auto_range = min,max,bset;
	self:update_range();
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:get_fields()
	if self.fields then
		return self.fields[ self.freq_idx ] ;
	end
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:get_shape(objs)
	if self.shape then
		return self.shape:copy();
	end
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:get_range()
	return self.range_min,self.range_max,self.auto_range;
end
--===========================================================================--

--===========================================================================--
--ui交互函数:刷新区间
function SBRSolver:refresh_range()
	if self.fields and self.fields[self.freq_idx] then  
		self.range_min,self.range_max = self.fields[ self.freq_idx ]:value_range();
		self:update_range();
	end
end
--===========================================================================--

--===========================================================================--
--ui交互函数:装载几何
function SBRSolver:load_geometry(path)
	print(path)
	Geometry.clear_all();
	--local g = Geometry.load_step(path);
	local g = Geometry.load(path);
	g:name('Loaded geometry')
	--Geometry.AddShape(g);
	g:build_shape();
	--g1 = g:surface_mesh();
	g:name('SBR')
	
	local bbox = g:bbox();
	bbox:extend( bbox:size()/10,bbox:size()/10,bbox:size()/10 )
	--print(bbox:xmin(),bbox:ymin(),bbox:zmin(),
					                               --bbox:xmax(),bbox:ymax(),bbox:zmax() )
	local b = Geometry.create_brick( {name='air box' ,pt1={bbox:xmin(),bbox:ymin(),bbox:zmin()}
										     ,pt2={ bbox:xmax(),bbox:ymax(),bbox:zmax()}
										     ,is_fill = false} );
	Geometry.AddShape(b);
	Geometry.AddShape(g)
	Gallop.View:perspective_view();
end
--===========================================================================--

--===========================================================================--
--ui交互函数:装载过滤几何
function SBRSolver:load_filter_geometry(path)
	print(path)
	--Geometry.clear_all();
	local obj = Geometry.object('SBRFilter');
	if obj then
		Geometry.RemoveShape(obj);
	end
	--local g = Geometry.load_step(path);
	local g = Geometry.load(path);
	g:name('SBRFilter')
	--Geometry.AddShape(g);
	g:build_shape();
	g:fill_color(1,0,0)
	--g1 = g:surface_mesh();
	--g1:name('SBRFilter')
	Geometry.AddShape(g)
	Gallop.View:perspective_view();
end
--===========================================================================--

--===========================================================================--
--ui交互函数:设置参数按钮
function SBRSolver:set_calculate_parameters()
	print( "set calculation parameters" );	
	self:create_calculation_model()
	self.binit = true;
	self.RCS = nil;
	--self.tranView:clearGraphs();
	self:plot();
	self:plot2D();
	self:update_progress_bar( 0 );
end
--===========================================================================--

--===========================================================================--
--创建SBR求解RCS模型
function SBRSolver:create_calculation_model()	

	Gallop.ParametersSweep:clear_all_results();
	Gallop.ParametersSweep:update_settings_from_ui();
	self:update_parameters_from_ui()
	self:init_parameters()		
	self:create_paras();
	self:set_paras();
	self:create_directions();

	---------------------------
	--设置透视视角
	Gallop.View:perspective_view();
	---------------------------		
	Geometry.synchronize();
end
--===========================================================================--

--===========================================================================--
--初始化参数
function SBRSolver:init_parameters()
		
	if self.polarHH==true then
		self.polars = {"HH"};
	end
	if self.polarVV==true then
		self.polars  = {"VV"};
	end
	if self.polarHH==true and self.polarVV==true then
		self.polars  = {"HH","VV"};
	end	

	self.probe_freq = nl.LinearSpace(self.minFreq,self.maxFreq,self.numFreq);
	if self.bFreqList==true then
		self.probe_freq = nl.new_vector(self.freqList);
	end		
	
	self.theta_angles = nl.LinearSpace(self.minTheta,self.maxTheta,self.numTheta);
	self.str_theta_angles = '('..self.minTheta .. ',' .. self.maxTheta .. ',' .. self.numTheta .. ')';
	if self.bThetaList==true then
		self.theta_angles = nl.new_vector(self.thetaList);
		self.str_theta_angles = self.thetaList;
	end
	
	self.phi_angles = nl.LinearSpace(self.minPhi,self.maxPhi,self.numPhi);
	self.str_phi_angles = '('..self.minPhi .. ',' .. self.maxPhi .. ',' .. self.numPhi .. ')';
	if self.bPhiList==true then
		self.phi_angles = nl.new_vector(self.phiList);
		self.str_phi_angles = self.phiList
	end
	
	self.polar_angles = nl.new_vector();
	for k,v in pairs(self.polars) do 
		if v == 'HH' then
			self.polar_angles:add(90);
		else
			self.polar_angles:add(0);
		end	
	end

	self.polarList = '';
	for k,v in pairs( self.polar_angles) do 
		self.polarList = self.polarList .. ' ' .. v ..';' ;
	end
end
--===========================================================================--

--===========================================================================--
--创建参数表
function SBRSolver:create_paras()
	Geometry.Parameters.add_variable("incident_angle","0","方位角: zx平面的夹角")
	Geometry.Parameters.add_variable("polar_angle","0",'极化角,垂直极化VV[0],水平极化HH[90]')
	Geometry.Parameters.add_variable("dfarbox","10",'远场盒子与几何结构距离')
	Geometry.Parameters.add_variable("dext","20",'入射源与几何结构距离')
	Geometry.Parameters.add_variable("elevation_angle","0",'俯仰角:与xy平面的夹角')	
end
--===========================================================================--


--===========================================================================--
--设置参数表
function SBRSolver:set_paras()

	local thetas = self.theta_angles;
	local phis = self.phi_angles;	
	local theta = thetas[1];
	local polar = self.polars[1];
	
	Geometry.Parameters.set_variable("polar_angle",tostring(self.polar_angles[1]))
	
	Geometry.Parameters.set_variable('incident_angle',''..self.theta_angles[1]);
	
	Geometry.Parameters.set_variable('elevation_angle',''..phis[1]);
	
	if self.bMonoRCS then 
		Geometry.Parameters.change_variable_sweep_state("incident_angle",true)
		Geometry.Parameters.change_variable_sweep_state("elevation_angle",true)
		Geometry.Parameters.change_variable_sweep_values("incident_angle",self.str_theta_angles)	
		Geometry.Parameters.change_variable_sweep_values("elevation_angle",self.str_phi_angles)
	else 
		Geometry.Parameters.change_variable_sweep_state("incident_angle",false)
		Geometry.Parameters.change_variable_sweep_state("elevation_angle",false)
	end

	Geometry.Parameters.change_variable_sweep_state("polar_angle",true)
	
	Geometry.Parameters.change_variable_sweep_values("polar_angle",self.polarList)

	Geometry.Parameters.valuate()
	
	--print( Gallop.ParametersSweep.index );
	
	Gallop.ParametersSweep:change_solver_index('SBRApp');
	
	Gallop.ParametersSweep:update_solver_list_to_ui();
	
	
		
	local preStr=[[	
-----------------------------------------------------------------------
--参数扫描前处理函数
function prerun(self)
	Gallop.SBRApp:initSolver();
end
-----------------------------------------------------------------------
]]

	local stepStr=[[	
-----------------------------------------------------------------------
--每一步扫描过程执行函数
function step(self,idx)
	--保存当前视图
	--self.view3D:saveSnapShot("view"..idx..".png")

	Gallop.SBRApp:preRunSolver(idx);

	--求解
	Gallop.SBRApp:runSolver();
	
	Gallop.SBRApp:postRunSolver(idx);	
		
end	
-----------------------------------------------------------------------
]]

	local postStr=[[	
-----------------------------------------------------------------------
--参数扫描后处理函数
function postrun(self)
	Gallop.SBRApp:collect_data();
	Gallop.SBRApp:update_border_to_ui();
	Gallop.SBRApp:update_freqlist_to_ui();
	Gallop.SBRApp:plot2D();	
	Gallop.SBRApp:export_data(Gallop.SBRApp.stemPath .. '-MonoRCS(SBR).dat' )
end
-----------------------------------------------------------------------
]]


	Gallop.ParametersSweep:set_solver_settings( stepStr  );
	Gallop.ParametersSweep:set_pre_settings( preStr);
	Gallop.ParametersSweep:set_post_settings( postStr);
	Gallop.ParametersSweep.data = nil;
		
	---------------------------
	--设置透视视角
	Gallop.View:perspective_view();
	---------------------------	
	
	Geometry.synchronize();

end
--===========================================================================--

--===========================================================================--
--初始化求解器
function SBRSolver:initSolver()
     self:get_path_from_project()
	
     collectgarbage()
     collectgarbage()
     self:create_sbr()     
     nl.num_threads(self.numThreads)

     print("size:", #self.theta_angles,#self.phi_angles,#self.probe_freq,#self.polars)
     self.RCS = nl.new_eigen_ctensor4(#self.theta_angles,
						       #self.phi_angles,
						       #self.probe_freq,
						       #self.polars);
     self.RCS:set_zero();
     
     if id_groups == 0 then 
		local path = unity.utf8_to_local(self.resultsPath .. "Parameters.txt")
		local file = io.open(path, "w")
		file:write('theta =\''..self.str_theta_angles..'\'\n');
		file:write('phi = \''..self.str_phi_angles ..'\'\n');
		file:write('polar = \''..self.polarList..'\'\n');
		file:close();
     end
end
--===========================================================================--

--===========================================================================--
--初始化求解器
function SBRSolver:preRunSolver(idx)
end
--===========================================================================--


--===========================================================================--
--求解后的后处理
function SBRSolver:runSolver(idx)
	collectgarbage()
	collectgarbage()
	self.sbr:build_shape()	
	if self.bMonoRCS then 
		self.tmp_RCS =  self.sbr:calculate_monoRCS_copolar_fields_mt( 1e9*self.probe_freq );
	else
		self.tmp_RCS =  self.sbr:calculate_bistaticRCS_RCS_mt( 1e9*self.probe_freq );
	end
end
--===========================================================================--

--===========================================================================--
--求解后的后处理
function SBRSolver:postRunSolver(idx)	
	if self.bMonoRCS then 
		self:postRunSolverMonoRCS(idx)
	else
		self:postRunSolverBiRCS(idx)	
	end
end

function SBRSolver:postRunSolverMonoRCS(idx)	
	print( "get mono RCS..." )	
	local index = Gallop.ParametersSweep.paraIndex[idx];
	local names = Gallop.ParametersSweep.names;
        local value = Gallop.ParametersSweep.paraList[idx]
	print('index:',index)
     	for kk=1,#self.probe_freq do
		self.RCS[{index[1]+1,index[3]+1,kk,index[2]+1}] =  self.tmp_RCS[kk]
	end	
	self:plot();
	local path = unity.utf8_to_local(self.resultsPath .. "MonoRCS-" .. idx ..".txt")
	local file = io.open(path, "w")
	for kk = 1,#names do 
		file:write(names[kk] .. '=' .. value[kk] .. '\n');
	end
	file:write('Freq[GHz]\tReal\tImag\tRCS(dBsm)\n')
	for kk=1,#self.probe_freq do
		file:write( self.probe_freq[kk] .. '\t' .. nl.real( self.tmp_RCS[kk] ) .. '\t' .. nl.imag( self. tmp_RCS[kk] ) .. '\t' .. self.data_post_fun(self.tmp_RCS[kk]) ..'\n');
	end
	file:close();
	self:update_progress_bar( idx / #Gallop.ParametersSweep.paraList * 100 );
end

function SBRSolver:postRunSolverBiRCS(idx)	
end
--===========================================================================--

--===========================================================================--
--初始化目录
function SBRSolver:get_path_from_project()
	self.projectPath = Geometry.get_project_path_utf8(  );
	assert( self.projectPath~='', "请先保存工程") 
	self.stemPath = Geometry.get_project_stem_utf8();
	
	self.path =  strippath(self.projectPath) .. '/';
	self.resultsPath = self.stemPath .. '-RCS/';	
	self.gapsPath = self.stemPath ..'-GAPs/';	
	self.name = getFilenameWithoutExt( self.projectPath );

	--self:print_paths();
	unity.mkdir(self.resultsPath)
	unity.mkdir(self.gapsPath)
end
--===========================================================================--

--===========================================================================--
--求解后的数据保存和读取
function SBRSolver:getResults(idx)	
	local data = nl.new_cvector(#self.probe_freq)
	--print(idx) 
	local index = Gallop.ParametersSweep.paraIndex[idx];
	for kk=1,#self.probe_freq do
		data[kk] = self.RCS[{index[1]+1,index[3]+1,kk,index[2]+1}] ;
	end	 
	return data;
end
function SBRSolver:setResults(data,idx)	
	local index = Gallop.ParametersSweep.paraIndex[idx];
	--print( #data, #self.probe_freq)
	for kk=1,#self.probe_freq do
		self.RCS[{index[1]+1,index[3]+1,kk,index[2]+1}] = data[kk];
	end	
end
--===========================================================================--

--===========================================================================--
--数据收集（主要针对MPI）
function SBRSolver:collect_data()	

	mpi.world.barrier();
	local reqs = mpi.request_list:new();
	local paraList = Gallop.ParametersSweep.paraList;
	local data = {};
	for k,v in pairs( paraList ) do
		if ( (k-1) % mpi.world.size() ) == mpi.world.rank()  then --如果是计算进程
			if mpi.world.rank() ~= 0 then	--如果不是主进程，则发送数据
				local index = Gallop.ParametersSweep.paraIndex[k];
				local rcs = nl.new_cvector(#self.probe_freq);
				for ii=1,#self.probe_freq do 
					rcs[ii] = self.RCS[{index[1]+1,index[3]+1,ii,index[2]+1}] ;
				end
				req = mpi.world.isend(0,k,rcs);
				reqs:push_back(req);				
			end
		else --如果不是计算进程则要收数据
			if mpi.world.rank() == 0 then	
				data[k] = nl.new_cvector();
				req = mpi.world.irecv( (k-1) % mpi.world.size(),k,data[k]);
				reqs:push_back(req);				
			end
		end
	end
	----------------------------------------------------------------------------------------------------------------------------------------
	--wait for all partical envaluation finish
	mpi.wait_all(reqs);
	mpi.world.barrier();	----------------------------------------------------------------------------------------------------------------------------------------
	
	if mpi.world.rank() == 0 then	
		for k,v in pairs(data) do 
			local index = Gallop.ParametersSweep.paraIndex[k];
			for ii=1,#self.probe_freq do 
				self.RCS[{index[1]+1,index[3]+1,ii,index[2]+1}]  = v[ii];
			end
		end	
	end
end
--===========================================================================--

--===========================================================================--
--创建入射角示意图
function SBRSolver:create_directions(bbox)
	self.bbox = Geometry.bbox();
	---------------------------
	--创建显示点
	local pts = Geometry.object("directions","shape");
	if pts == nil then 
		pts = Geometry.NewPointSet("directions");
		Geometry.AddShape(pts);
	end
	pts:clear();
	local bbox = self.bbox;
	--local pts = Geometry.NewPointSet();
	local s = bbox:size()*0.5;
	local xc = (bbox:xmin()+bbox:xmax())/2.0
	local yc = (bbox:ymin()+bbox:ymax())/2.0
	local zc = (bbox:zmin()+bbox:zmax())/2.0
	
	eang = 0;
	for _,phi in pairs (self.phi_angles) do 
		local eang = phi * pi/180
		for _,ang in pairs(self.theta_angles) do		
			pts:insert( s * math.cos( ang *pi/ 180 ) * math.cos( eang )+ xc,
					s * math.sin( ang *pi/ 180 ) * math.cos( eang )+ yc,
					s * math.sin( eang ) + zc,
				       -s/20*math.cos( ang *pi/ 180 )  * math.cos( eang ), 
				       -s/20*math.sin( ang *pi/ 180 )  * math.cos( eang ),
				       -s/20 * math.sin( eang )	);
		end
	end
	pts:point_size(8)
	pts:fill_color(1,0,0)
	
	local axis = Geometry.object("incident","shape");
	if axis == nil then 
		axis = Geometry.NewAxis({name="incident",xName="E",yName="H",zName="k"});
		Geometry.AddShape(axis);
	end
	axis:load_identity_matrix();
	axis:scale(s/5,s/5,s/5)
	axis:rotate(zAxis,'polar_angle');
	axis:rotate(yAxis,-90);
	axis:translate(s,0,0);
	axis:rotate(yAxis,'-elevation_angle');	
	axis:rotate(zAxis,'incident_angle');	
	axis:translate(xc,yc,zc);
end
--===========================================================================--

--===========================================================================--
--初始化参数
function SBRSolver:create_sbr()
	
	self.shape = Geometry.object('SBR');	
	assert(self.shape, '未找到命名为SBR的物体')
	
	--新建一个射线追踪器
	if self.sbr == nil then
		self.sbr = Geometry.object('SBRSolver');
		if self.sbr == nil then 
			self.sbr = Geometry.NewSBR('SBRSolver')		
			self.sbr:is_outline(false)
			self.sbr:is_point(false)
		end 
	end
	self.sbr:mesh(self.shape:copy())
	if  Geometry.object('SBRSolver') == nil then 
		Geometry.AddShape(self.sbr);
	end
	
	local vec = Geometry.NewVector()
        vec:xyz("-1*sin((90-elevation_angle)*pi/180)*cos(incident_angle*pi/180)","-1*sin((90-elevation_angle)*pi/180)*sin(incident_angle*pi/180)","-1*cos((90-elevation_angle)*pi/180)")         
        self.sbr:incident_kvector(vec);
        local vec1 = Geometry.NewVector()

	vec1:xyz("cos((90-elevation_angle)*pi/180)*cos(incident_angle*pi/180)*cos(polar_angle*pi/180)-sin(incident_angle*pi/180)*sin(polar_angle*pi/180)",
	"cos((90-elevation_angle)*pi/180)*sin(incident_angle*pi/180)*cos(polar_angle*pi/180)+cos(incident_angle*pi/180)*sin(polar_angle*pi/180)",
	"-1*cos(polar_angle*pi/180)*sin((90-elevation_angle)*pi/180)")
	self.sbr:incident_Hvector(vec1)
       --------------HH极化
        --vec1:xyz("-1*sin(incident_angle*pi/180)","cos(incident_angle*pi/180)","0")
	--self.sbr:incident_Hvector(vec1)
	--设置射线管密度
	self.sbr:distance_of_neighboring_rays(self.minMesh)
	self.sbr:reserve_traces(self.bReserveTraces)
	self.sbr:has_filter(self.bFilter)
	self.sbr:is_inner_filter(self.bInnerFilter)
	if self.bFilter==true then
		local shape1 = Geometry.object('SBRFilter');	
		if shape1 then
			--print('set filter_obj' )
			self.sbr:filter_obj(shape1:copy());
		end
	end	
end
--===========================================================================--

--===========================================================================--
--初始化数据
function SBRSolver:init_data()
	self.data = nl.new_eigen_matrix(#self.theta_angles,#self.probe_freq+2)		
	for ii=1,#self.theta_angles do
		for jj=1,#self.phi_angles do
		end
	end
end
--===========================================================================--

--===========================================================================--
--计算
function SBRSolver:run_pond()
	Gallop.ParametersSweep:run_sweep();
end
--===========================================================================--

--===========================================================================--
--实时更新1D图片
function SBRSolver:plot()
	if Gallop.console==true then return end;
	self.rcsView:clearGraphs();
	if self.RCS == nil then 
		self.rcsView:rescaleAxes();
		return; 
	end;
	for p, pp in pairs (self.polars ) do 
		for k,v in pairs( self.probe_freq ) do 
			val = nl.new_vector(#self.theta_angles);
			for jj=1,#self.phi_angles do 
				for ii = 1,#self.theta_angles do 
					val[ii] = self.data_post_fun( self.RCS[{ii,jj,k,p}] )
				end
				self.rcsView:addGraph( self.theta_angles, val, self.probe_freq[k] ..' GHz 俯仰角 '..self.phi_angles[jj] ..'° ' .. pp)
			end
		end 
	end
	Gallop.View:set_plot_lable( self.rcsView, '方位角 [°]', self.ylable_names[self.scale_type] );
	self.rcsView:rescaleAxes()
end

function SBRSolver:add_theta_graph(phi,freq,polar)
	local val = nl.new_vector(#self.theta_angles);
	for ii = 1,#self.theta_angles do 
		val[ii] = self.data_post_fun( self.RCS[{ii,phi,freq,polar}] )
	end
	self.rcsView:addGraph( self.theta_angles, val, self.probe_freq[freq] ..' GHz 俯仰角 '..self.phi_angles[phi] ..'° ' .. self.polars[polar])
end
--===========================================================================--

--===========================================================================--
--实时更新2D云图图片
function SBRSolver:plot2D()
	if Gallop.console then return end;
	self.fields = {};
	if self.RCS == nil then 
		self:update_view();
		return; 
	end;
	if self.shape== nil then 
		self.shape = Geometry.object('SBR');
	end	
	if self.shape == nil then
		self:update_view();
		return; 
	end;

	local idx = 0;
	self.bbox = self.shape:bbox();
	for ifreq = 1,#self.probe_freq do 
		for ipolar = 1,#self.polar_angles do 
			idx = idx+1;
			local fields = self:create_fields();
			for itheta = 1,#self.theta_angles do 
				for iphi = 1,#self.phi_angles do 
					local f = nl.new_vector(1);
					f[1] =  self.data_post_fun( self.RCS[{itheta,iphi,ifreq,ipolar}] )
					fields:insert((90-self.phi_angles[iphi])/180*pi,self.theta_angles[itheta]/180*pi,f)
				end
			end
			
			fields:scale(self.bbox:size()*1)
			fields:translate(self.bbox:xcenter(),self.bbox:ycenter(),self.bbox:zcenter())
			fields:set_range(-50,0)
			cm = fields:colormap()
			cm:min(-50);
			cm:max(0);
			fields:build_polyhedron();
			fields:build_shape();
			self.fields[idx] =fields;
		end
	end
	--self:change_selected_freq(0)
	--self:update_view();
	self:update_shape();
	self:init_view2D();
	
	self.freq_idx =1;
	self:update_freqlist_to_ui();
	self:update_freq_index_to_ui();
	self:update_fields_to_view();
	self:update_range()
	self:update_view();
end
--===========================================================================--

--===========================================================================--
--格式化数据
function SBRSolver:get_rcs_data()
	local RCS = self.RCS;
	local nphi = #self.phi_angles;
	local ntheta = #self.theta_angles;
	local npolar = #self.polar_angles;
	
	if RCS~=nil then 
		local header = nl.new_svector();
		header:add('Index')
		header:add('Theta(°)')
		header:add('Phi(°)')
		header:add('Polar angle(°)')
		for ifreq = 1,#self.probe_freq do  			
			header:add( 'real(Farfields)@'..self.probe_freq[ifreq]..' GHz' )
			header:add( 'imag(Farfields)@'..self.probe_freq[ifreq]..' GHz' )
		end
		local data = nl.new_eigen_smatrix(#self.theta_angles*#self.phi_angles*#self.polar_angles+1,2*#self.probe_freq+4)
		data[{1,1}] = 'Frequency[GHz]';
		data[{1,2}] = '--';
		data[{1,3}] = '--';		
		data[{1,4}] = '--';		
		for ifreq = 1,#self.probe_freq do  
			data[{1,ifreq+4}] = tostring( self.probe_freq[ifreq] );
		end
		idx = 0;
		for ipolar = 1,#self.polar_angles do 
			for iphi = 1,#self.phi_angles do 
				for itheta = 1,#self.theta_angles do 
					idx = idx+1;
					--local idx = ntheta * (iphi-1) + itheta+1;
					data[{idx,1}] = tostring( idx-1 );
					data[{idx,2}] = tostring( self.theta_angles[itheta] );
					data[{idx,3}] = tostring( self.phi_angles[iphi] );
					data[{idx,4}] = tostring( self.polar_angles[ipolar] );
					for ifreq = 1,#self.probe_freq do  
						data[{idx,ifreq*2+3}] = tostring( nl.real( RCS[{itheta,iphi,ifreq,ipolar}] ) );
						data[{idx,ifreq*2+4}] = tostring( nl.imag( RCS[{itheta,iphi,ifreq,ipolar}] ) );
					end
				end
			end
		end
		return data,header;
	end 
end
--===========================================================================--

--===========================================================================--
--格式化数据
function SBRSolver:get_rcs_data()
	local RCS = self.RCS;
	local nphi = #self.phi_angles;
	local ntheta = #self.theta_angles;
	local npolar = #self.polar_angles;
	
	if RCS~=nil then 
		local header = nl.new_svector();
		header:add('Index')
		header:add('Theta(°)')
		header:add('Phi(°)')
		header:add('Polar angle(°)')
		for ifreq = 1,#self.probe_freq do  		
			if self.scale_type ==3 or self.scale_type==4 then
				header:add( 'real(Farfields)@'..self.probe_freq[ifreq]..' GHz' )
				header:add( 'imag(Farfields)@'..self.probe_freq[ifreq]..' GHz' )
			end
			if self.scale_type == 1 then
				header:add( 'RCS(dBsm)@'..self.probe_freq[ifreq]..' GHz' )
			end
			if self.scale_type == 2 then
				header:add( 'RCS(sm)@'..self.probe_freq[ifreq]..' GHz' )
			end
		end
		local data;
		if self.scale_type ==3 or self.scale_type==4 then		
			data = nl.new_eigen_smatrix(#self.theta_angles*#self.phi_angles*#self.polar_angles+1,2*#self.probe_freq+4)
		else
			data = nl.new_eigen_smatrix(#self.theta_angles*#self.phi_angles*#self.polar_angles+1,#self.probe_freq+4)
		end
		data[{1,1}] = 'Frequency[GHz]';
		data[{1,2}] = '--';
		data[{1,3}] = '--';		
		data[{1,4}] = '--';		
		for ifreq = 1,#self.probe_freq do  
			data[{1,ifreq+4}] = tostring( self.probe_freq[ifreq] );
		end
		idx = 0;
		for ipolar = 1,#self.polar_angles do 
			for iphi = 1,#self.phi_angles do 
				for itheta = 1,#self.theta_angles do 
					idx = idx+1;
					--local idx = ntheta * (iphi-1) + itheta+1;
					data[{idx,1}] = tostring( idx-1 );
					data[{idx,2}] = tostring( self.theta_angles[itheta] );
					data[{idx,3}] = tostring( self.phi_angles[iphi] );
					data[{idx,4}] = tostring( self.polar_angles[ipolar] );
					for ifreq = 1,#self.probe_freq do
						if self.scale_type ==3 or self.scale_type==4 then		
							data[{idx,ifreq*2+3}] = tostring( nl.real( RCS[{itheta,iphi,ifreq,ipolar}] ) );
							data[{idx,ifreq*2+4}] = tostring( nl.imag( RCS[{itheta,iphi,ifreq,ipolar}] ) );
						end
						if self.scale_type == 1 or self.scale_type== 2 then
							data[{idx,ifreq+4}] = tostring( self.data_post_fun( RCS[{itheta,iphi,ifreq,ipolar}] ) );
						end
					end
				end
			end
		end
		return data,header;
	end 
end
--===========================================================================--

--===========================================================================--
--查看数据
function SBRSolver:view_rcs_data()
	local data,header = self:get_rcs_data()
	if data and header then
		self.rcsdataView  = Gallop.DataPost:create_data_table_to_ui();
		self.rcsdataView:set_data(data);
		self.rcsdataView:set_header(header);
	end 
end
--===========================================================================--

--===========================================================================--
--导出数据
function SBRSolver:export_data(utf8path)
	--self.rcsView:exportData(path)
	local localpath = unity.utf8_to_local(utf8path);
	local data,header = self:get_rcs_data()
	local file = io.open(localpath,'w');
	if data and header then
		for k,v in pairs(header) do 
			file:write( v .. '\t' );
		end
		file:write( '\n' )
		local nrows = data:rows();
		local ncols = data:cols();
		for ii=1,nrows do 
			for jj=1,ncols do 
				file:write( data[{ii,jj}] .. '\t' );
			end
			file:write( '\n' )
		end
	end
	file:close();
end
--===========================================================================--

--===========================================================================--
--画单个状态的RCS图
function SBRSolver:plot_theta_graph(phi,freq,polar)
	if Gallop.console==true then return end;
	self.rcsView:clearGraphs();
	if self.RCS == nil then 
		self.rcsView:rescaleAxes();
		return; 
	end;
	self:add_theta_graph(phi,freq,polar);
	Gallop.View:set_plot_lable( self.rcsView, '方位角 [°]', self.ylable_names[self.scale_type]);
	self.rcsView:rescaleAxes()
end
--===========================================================================--

--===========================================================================--
--导出报告所用的数据
function SBRSolver:save_reports()
     self:get_path_from_project()
     self.reportsPath = self.stemPath ..'-Reports/';	
     unity.remove(self.reportsPath);
     unity.mkdir(self.reportsPath);
     print(self.scale_type)
     self:change_scale_type(2)
	for p, pp in pairs (self.polars ) do 
		for k,v in pairs( self.probe_freq ) do 
			for jj=1,#self.phi_angles do 
				self:plot_theta_graph(jj,k,p);
				--self.rcsView:apply_function('透过率频谱');
				self.rcsView:apply_function('General RCS');
				pathstr = self.reportsPath .. '频率' .. self.probe_freq[k] ..' GHz_俯仰角 '..self.phi_angles[jj] ..'°_' .. self.polars[p].. '.png';
				self.rcsView:savePng( unity.utf8_to_local(pathstr) )
			end
		end 
	end
	
	for kk=1,#self.probe_freq * #self.polar_angles do 
		self:change_selected_freq(kk)
		self.view3D:saveSnapShot( self.reportsPath .. '/' .. kk .. '.png' )
	end
     self:cal_generic_rcs()
     self:save_word()

end
--===========================================================================--


--==============================================================================================================--
function cal_generic_rcs(x1,y1)
	
	local sum = 0;
	local num = 0;
	for ii=1,#x1 do
		if nl.abs( nl.wrap_phase_deg(x1[ii]*4)/4 ) > 4.5 then
			sum = sum+y1[ii]
			num = num + 1;
		end
	end
	sum = sum/num
	return sum;
end
--==============================================================================================================--

--==============================================================================================================--
function SBRSolver:add_style(doc,name)
	local newStyle = doc.Styles:Add(name ,1) -- 1 代表段落样式
	newStyle.BaseStyle = "正文" -- 基于“Normal”样式
	newStyle.Font.Name = "黑体" -- 设置字体
	newStyle.Font.Size = 12 -- 设置字号
	newStyle.ParagraphFormat.Alignment = 1 -- 设置对齐方式，1 代表居中
end
--==============================================================================================================--

--==============================================================================================================--
function SBRSolver:add_fig_to_doc(doc,imagePath,label)

	insertionPoint = doc.ActiveWindow.Selection
	--insertionPoint.Style = "Heading 1"
	insertionPoint.Font.Name = "Tahoma"
	local function replaceBackslash(inputStr)
	    local outputStr = string.gsub(inputStr,  "/","\\")
	    return outputStr
	end

	insertionPoint:TypeParagraph()
	insertionPoint.Style = "图"
	insertionPoint.InlineShapes:AddPicture(replaceBackslash(imagePath))

	insertionPoint:TypeParagraph()
	insertionPoint.Style = "图示"
	insertionPoint:TypeText( label )
end
--==============================================================================================================--

--==============================================================================================================--
function SBRSolver:add_all_fig_to_doc(doc)
	for p, pp in pairs (self.polars ) do 
		for k,v in pairs( self.probe_freq ) do 
			for jj=1,#self.phi_angles do 
				pathstr = self.reportsPath .. '频率' .. self.probe_freq[k] ..' GHz_俯仰角 '..self.phi_angles[jj] ..'°_' .. self.polars[p].. '.png';
				label = 'RCS角度谱( @频率' .. self.probe_freq[k] ..' GHz_俯仰角 '..self.phi_angles[jj] ..'°_' .. self.polars[p].. ')';
				print( label )
				self:add_fig_to_doc( doc, pathstr, label );
			end
		end 
	end
end
--==============================================================================================================--

--==============================================================================================================--
function SBRSolver:cal_generic_rcs()
	local w = Gallop.DataPost:create_data_table_to_ui();
	local data = nl.new_eigen_smatrix(#self.probe_freq*#self.phi_angles*#self.polar_angles,4)
	idx = 1;
	self:update_post_fun();
	for p, pp in pairs (self.polars ) do 
		for k,v in pairs( self.probe_freq ) do 
			for jj=1,#self.phi_angles do 
				local val = nl.new_vector(#self.theta_angles);
				for ii = 1,#self.theta_angles do 
					val[ii] = self.data_post_fun( self.RCS[{ii,jj,k,p}] )
				end
				rcs = cal_generic_rcs( self.theta_angles, val );
				data[ {idx, 1} ] = tostring(pp);
				data[ {idx, 2} ] = tostring(v);
				data[ {idx, 3} ] = tostring(self.phi_angles[jj]);
				data[ {idx, 4} ] = tostring(rcs);
				idx = idx + 1;
			end
		end 
	end
	w:set_data( data );
	self.generic_rcs_data  = data;
	local header = nl.new_svector();
	header:add('Polar')
	header:add('频率(GHz)')
	header:add('仰角(°)')
	header:add('Generic RCS(sm)') 
	self.generic_rcs_header = header
	w:set_header( header );
end
--==============================================================================================================--

function SBRSolver:get_rcs_elm_data(r,c)
	return tostring(self.generic_rcs_data[{r,c}])
	--return 'adf'
end

--==============================================================================================================--
function SBRSolver:save_word()
	self:generate_doc(self.reportsPath..'/a.docx');
	require("luacom")
	-- 启动Word应用程序
	local word = luacom.CreateObject("Word.Application")
	word.Visible = true -- 让Word应用程序可见
	-- 添加一个新文档
	local doc = word.Documents:Add()
	self:add_style(doc,'图')
	self:add_style(doc,'图示')
	-- 获取文档的范围
	local content = doc.Content	-- 插入表格，参数分别是行数和列数
	
	-- 插入表格，参数分别是行数和列数
	local table = doc.Tables:Add(content, self.generic_rcs_data:rows()+1, self.generic_rcs_data:cols() ) -- 插入一个3x3的表格
	-- 填充表格数据
	for r = 1, self.generic_rcs_data:rows() do
	    for c = 1,  self.generic_rcs_data:cols() do
		local cell = table:Cell(r+1, c)
		cell.Range.Text =  self.generic_rcs_data[{r,c}]
	    end
	end
	
	local cell = table:Cell(1, c)
        for c = 1,  self.generic_rcs_data:cols() do
	    local cell = table:Cell(1, c)
	    cell.Range.Text =  self.generic_rcs_header[c]
       end
	
	-- 应用样式到表格。确保样式名称与 Word 中的实际样式名称匹配。
	table.Style = "网格型"
	
	-- 将选择区域移动到表格后面
	selection = doc.ActiveWindow.Selection
	selection:SetRange(table.Range.End + 1, table.Range.End + 1)
	selection:Collapse(0) -- 移动到区域的结束位置
	
	self:add_all_fig_to_doc(doc);

	-- 保存文档
	doc:SaveAs(self.reportsPath..'rcs.docx')
	-- 关闭文档
	doc:Close(false) -- 参数false表示关闭文档时不保存更改，如果需要保存更改，请传递true
	 --退出Word应用程序
	word:Quit()
	 --释放COM对象
	word = nil
end
function SBRSolver:save_word()
     collectgarbage()
     collectgarbage()
          collectgarbage()
     collectgarbage()
	self:generate_doc(self.reportsPath..'/a.docx');
end
--==============================================================================================================--

return SBRSolver;